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FOREWORD and SUMMARY 


This document is the final report to the National Aeronautics and Space 
Administration, George C. Marshall Space Flight Center, Marshall Space Flight 
Center, Huntsville , Alabama. This effort was performed by Continuum, Inc., 
Huntsville, Alabama under Contract Number NAS 8-349^6. The technical 
representative from MSFC was Mr. Klaus Gross. 

The information presented here deals with a different approach to modelling of 
•the thermochemistry of rocket engine combustion phenomena. While the approach 
is essentially unverified by extensive test data, it does advance a rationale 
that explains some troubling discepancies in the conventional model. The 
methodology described here is based on the Continuum, Inc. hypothesis of a new 
variational principle applicable to compressible fluid mechanics. This 
hypothesis is extended to treat the thermochemical behavior of a reacting 
(equilibrium) gas in an open system. 



NOMENCLATURE 




The nomenclature used in this report is consistent with that given in 
reference (2) with the following additions; 

♦qj forcing function modifying the Gibbs function 

E 

^ defined in text 

defined in text 


. n 


defined in text 
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This study deals with the determination of chemical equilibrium in steady 
flowing systems subject to downstream boundary conditions. Research at 
Continuum, Inc. has recently lead to the hypothesis of a new variational 
principle for compressible fluid mechanics. This principle , which has 
permitted impressive advances in numerical computations of unsteady flows, is 
extended to treat the one dimensional combustion and expansion of reacting 
gases in rocket engines. 

. The study of the flow of reacting gases has perhaps its most important 
application in rocket engines for use in the aerospace field. The performance 
of an engine in these applications is a critical consideration. For this 
reason precise methods are needed to predict, in a preliminary design, the 
thrust and mass flow of a proposed engine. Unfortunately this performance 
evaluation is extremely complex. Many imperfectly understood phenomena combine 
to create the combustion process while the ensuing hostile environment makes 
experimental measurements of all but gross quantities virtually impossible. 

The concept of equilibrium chemical reactions is invaluable in the performance 
evaluation of rocket engines. Although all reactions in a real engine progress 
at a finite rate, the assumption of equilibrium (infinite rates) provides 
meaningful information with a comparatively simple calculation. The 
equilibrium methodology is well acccepted but, because of its importance, a 
reexamination of the concepts is justifiable. 

Another reason for concern is that, after all effects such as droplet 
vaporization, reaction kinetics, fuel striations, divergence losses, boundary 
layer losses, etc. are accounted for, unacceptable errors may still exist in 
the performance evaluation. This discrepancy has been termed a combustion 
efficiency problem. The explanation usually given is that some of the above 
mechanisms may be imperfectly understood and that not enough experimental data 
has been gathered to adequately define some of the required parameters. 
Although this argument has merit and is certainly correct to some degree, 
other phenomena which would account for the errors are possible. 
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Several factors influenced the decision to pursue the questioning of 
equilibrium concepts. Performance calculations based conventional equilibrium 
concepts seemed to give systematically erroneous results. The adiabatic flame 
and isentropic expansion solution, typified by Zeleznick and Gordon gave 
good mass flows but overpredicted thrust and, therefore, specific Impulse. 


Now the adiabatic flame calculation assumes an infinite area duct at the 
combustion front. Real engines, of course, posses a finite area duct and it is 
well known that combustion in such a duct causes a total pressure drop. This 
effect is discussed in many texts using simple heating relations. 

.Consideration of this effect in rocket performance analyses served to reduce 
the thrust to approximately the right value but had the disconcerting effect 
of reducing the mass flow rate proportionately such that the impulse remained 
in error. 


The chemical equilibrium calculation itself had never been questioned. 
Proceeding from the premise that the effects were related and the apparent 
discrepancy of mass, thrust and impulse showed up in even the most basic 
building block (the equilibrium analysis) of the performance evaluation, it 
was postulated by Prozan that the throat choking condition acts as a 

constraint to the combur.tion process. That study was not well received and the 
matter remained unquestioned and unresolved for many years. 

Recently, however, research at Continuum led to the hypothesis of a 
variational principle for compressible flow. Prozan advanced the 

following statements for •consideration as the long sought after variational 
principle governing fluid mechanics; 

I. For the maximum stability of any real system, which may or may not be 
isolated from a larger universe, the rate of production of extensive entropy 
is a maximum subject to the constraints upon its* formation. 

II. For the equilibrium or time invariant state of such a system, the 
extensive entropy production must be a maximum subject to the constraints upon 
its formation. Consequently the entropy delivered to the universe, should the 
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system not be isolated, must also be a maximum. 




The first statement has been successfully used in unsteady flow of a 
compressible ideal gas to produce a numerical analog which is absolutely 
conservative and stable. The hypothesis, then, has been supported in at least 
one branch of fluid mechanics. Since the original study of reference is 
based on the same general principle, it seemed possible, even probable, that 
the thermochemical behavior of flowing systems was related to the hypothesized 
principle. Now the second statement is a consequence of the first statement 
for steady state. It states that the delivery of extensive entropy to the 
universe is a maximum for an open system. 


It is useful at this juncture to review the premise upon which conventional 
thermochemistry calculation arc baaed. The aLatementa of Gibba ' provide the 
foundation for the analysis of reference (1) and similiar calculations. The 
classical Gibbs statements are; 


I. "For the equilibrium of any isolated system it is necessary and sufficient 
that in all possible variations of the state which do not alter its energy, 
the variation of its entropy shall vanish or be negative." 

II. "For the equilibrium of any isolated system it is necesssary and 
sufficient that in all possible vaiations in the state of the system which do 
not alter its entropy, the variation of its energy shall either vanish or be 
negative." 

The author has no argument with the Gibbs statements. The point of contention 
is the interpretation of these statements and their application to 
thermochemical calculations in open systems. First notice that Gibbs 
statements deal with Isolated systems. Gibbs therefore dealt with 
thermostatics but not necessarily with thermodynamics. Secondly the Gibbs 
statements do refer to all possible variations in the state. It is the 
contention of this study that any final solution value which does not satisfy 
the throat condition is not a possible state. 

A final point which should be mentioned before proceeding with the discussion 
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relates to the formation of the Gibbs funotional. The previous two statements 
by Gibbs and the assumption of nonstant temperature and pressure permit one to 
develop the point thermodynamic function known as Gibbs free energy. It is 
this free energy function which is used in conventional thermodynamic 
calculations. The consequences of considering constant temperature and 
pressure in a changing pressure, temperature situation need to be addressed. 

The following technical discussion closely follows the previous study found in 
reference (2). The equilibrium combustion constrained by an ideal gas 
expansion to the throat is duplicated from that reference. In addition an 
equilibrium expansion from chamber to throat is included. In the interest of 
.brevity the salient features of that reference will not be repeated, but 
alterations and expansions upon that study will be discussed. Whereever 
possible, the same nomenclatijre and units are used. 
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TECHNICAL DISCUSSION 


Ideal Expansion from Chamoer to Throat 


The procedure found ui reference (?) is followed exactly in this study with 

the exception that the first two equations found in equation set (23) of that 

M T 

report were solved analytically rather than numerically for and , This 
led to a final set of equations to solve which are much simpler than those 
shown in the reference. This new set follows; 

n 


p ♦ b^v 


P 


inj 


ba 2 

h — V : ho 
2 


In (A*/A^) = 0 


where 


E 

I (h /T - il; — ) 
3T 


and where 






1 


Note that the system of equations has been cast so that it looks like the 
normal Gibbs function treatment except for the throat constraint and the new 
forcing function . The results were in substantial agreement with the 
reference (differing slightly because different thermodynamic parameters are 
available today). 
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In the equilibriam expansion from the chamber to the throat, the same 
optimization function was used. That is, the entropy was maximized subject to 
the system constraints. In this treatment, however, instead of varying a 
Lagrange multiplier to satisfy the throat constraint, the number of moles of 
the species OH was used to satify the throat constraint. This is actually the 
amount of OH at the throat rather than the amount of post combustion OH. 
Proceeding in the same fashion as the reference the entropy maximization 
statement was constrained by momentum and energy and the throat choking 
condition. The other difference between the equilibrium expansion and the 
.ideal gas expansion comes from the integration of the equations of motion 
during the expansion rather than the previous ideal treatment. For each 
station during the expansion or integration the entropy is maximized subject 
to the momentum and energy constraints locally. The system of equations used 
to describe the process is the same as those shown earlier except that; 


p 2 3h 2 

r = ( bjv fti-C2TJl2/4))/(T — *-b2V Qi) 

k-1 k 

(a‘‘ + A*^) 

whereJli = pn/(pn-biv) and ^ 2 = b,v/(prrhiv) where q = r — 

2 A*^ 


W -.1 

In the solution procedure the chamber solution is made by setting A"' and A 
(denoting area of the previous station and the current station respectively) 
to the combustor area. The number of moles of OH is held fixed. The 
Integration is then performed from chamber to tnroat by decreasing the 
pressure slightly from station to station. The crossectional area is found 
from the momentum equation for Lms new pressure. 


It was very interesting to note that the expansion was isentropic even though 
the entropy was being maximized. This was what was expected, but attaining an 
isentropic expansion with this procedure was never the less encouraging. The 
seemingly contradictory result can be interpreted in the following fashion. 

The variational principle states that the entropy delivered to the universe 
must be a maximum. This must be the case even during an isentropic expansion. 
The locus of states that represents equilibrium is a series of maximum values 
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of entropy whoso value during; the expansion is constant. This concept is not 
really hard to understand since, if the entropy is truly a maximum, then 
during a slow and continuous expansion it should stay const. iut or isentropic. 
The principle, if correct, should predict this behavior. No success has 
resulted from a limited attcfmpt to apply the principle beyond the throat. To 
solve this problem was beyond the scope of work of the contract and beyond the 
available funds in the contract. 


A computer program has bf*cn writt#?n to perform the calculations of equilibrium 
combustion and equilibrium orideal expansion between combustor and throat. 

This program is included in this report as an Appendix. Hopefully the 
.calculations are sufficiently well anotated that the program is easily 
understood . 

Where ever possible common terminology is used between report and program. 


Conclusions 


The entire subject is of great importance and certainly Is of great academic 
interest. The corusipl.s appear t,o h.ive potential to explain some* of the 
propulsion prediction discrepancies but are essentially unverified. The issue 
is best resolved by a dedicated test program which attempts to experimentally 
provide data for cases in which the constrained theory will predict a large 
difference from the conventional approach. 

Before such a test program be conducted, how»«ver, the analysis should be 
extended into th« supersonic regime of the nozzle. A successful solution of 
supersonic expansion would then allow one to more accurately assess the 
performance alterations caused by the new modelling concept. 
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ooo 


100 CALL HASA 
ICNST-7 

C COMPUTE FORClNCi FUNCTION fO ALTER GIBBS FUNCTION 
DO 110 1*1,6 

PHI<I>*ALAME*<HA< I )-I»W*L»HDr> 






110 CONTINUE 

C FIND NEW EQUILIBRIUM POSITION FOR THIS P.T 
CALL EQUILIICNST) 

C SOLVE MOMENTUM ANB ENERGY FOUATIONS FOR NEW P»T 
CALL HASA 
BM«CS#C1/(AC*AC:> 

BE*C4*C1*C 1 / ( AC«AC > 

P». 5*PF+SQRT ( . 25#PF»PF-BM*T/ W ) 

ENU-T/(P*W> 

HERR«H4- . 5*BE*ENU#ENtl-HT 
C DRIVE ENERGY ERROR TO /rR‘0 

•CALL ITSUB<T,HERR,20,. l.SAVEIT) 

IBR“IFIX<SAVEIT( 1 ) > 

GOTO i 100, K»0, 100, K»0, .' uO, 20«.0 , IBR 
200 CONTINUE 

C SOLUTION OF ALL EQUATIONS EXftPT THE THROAT CONSTRAINT HAS 
C CONVERGED FOR TH»S VAUJF iTT- At .AMF 
EM-SQRT < BM»T / < W*UAMMA ) > / P 
TFOT-1 . 4 . 5iM GAMMA-1 .) *EM#FM 
C FIND THE CHOKING AREA FOR THIS ALAME 
ASOACL-ATDET (GAMMA, EM. AsTAR > 

AS0ATL«AS0ACL4AL0G(AC/AT ) 

C CHANGE ALAME TO DRIVE L00< ASTAR/AT ) TO ZERO 
CALL ITSUBC ALAME, ASOATL . ?0, . 0001 , SAVE II > 

IBR«IFIX(SAVEIL(1>) 

GOTO (20,20,20.20,300. 300). IBR 
300 CONTINUE 

C SOLUTION HAS CONVERGED - PRINT COMBUSTOR AND THROAT RESULTS 
RH0«P*W/(C1«T) 

VEL*1./(RHCI*AC) 

WRITE (6. 9930) 

URITE(6*9800> 

HRtTE(S.9810> T.P.RHO.EM. VEL 
WR1TE(S,9B40) W 
WRITE(6,9410> 

WR1TE(6,9820> (XA( I > . 1*1 ,4) 

MR1TE(S*9830> (EMUA(l).I 1,6) 

WR1,TE(S*9100) 

TS-T*TFOT/( .5# (GAMMAS 1 . > ) 

P8-P*f ITS/T >**<0AMMA/ (GAMMA- 1 .)> ) 

RH0S>iP8»W/(Cl*TS) 

VCL8»1./(I»I0S*AT> 

A18P8«IP8*AT-h:3*VELS) / (9. 307#C3# 

MatEl4f9900> 

UR1TCC6»9910> AISPS.TS.PS.RHOS.VELS 
URlTK(4tPI00> 


ENTBI TIC CALCULATION OF THf EQUILIBRIUM EXPANSION TO THE THROAT 



ClflO 


100 CALL NASA 
ICNST»7 

C C0W>UTE FORCING FUNCTION TO ALTER GIBBS FUNCTION 
DO 110 1*1,6 

PHI ( I ) «ALAME* ( HA ( I ) -T frW*UHD X ) 

110 CONTINUE 

C FIND NEW EQUILIBRIUM POSITION FOR THIS P,T 

CALL EQUIL< ICNST) 

C SOLVE MOMENTUM AND ENERGY FOLIATIONS FOR NEW P,T 
CALL HASA 
BM=C3*C1 / ( AC* AC ) 

BE*C4*C1 *C 1 / ( AC*AC ) 

P*. 5*PF+SQRT ( . 25*F'F*PF-Bf1*T/W ) 

ENU*T/<P*W) 



FCRR»H+. 5*BE*ENU*ENLI--HF 
C DRIVE ENERGY ERROR TO 7PR0 

tW-L I TSUB < T , HERR , 20 , . 1 , SAVE I T ) 

IBR=IFIX(SAVFIT( 1 ) ) 

GOTO <100,100,1 00 , 1 00 , IM* ) , 200 ) , I BR 
200 CONTINUE 

C SOLUTION OF ALL EQUATIONS EXCEPT THE THROAT CONSTRAINT HAS 
C CONVERGED FOR THIS VALUF np Al. AME 
EM=SQRT < BM*T / ( W* GAMMA ) .» / 1 ' 

TFOT* 1 . + . 5* ( GAMMA - 1 . ) *EM* FM 
C FIND THE CHOKING AREA FOR THIS ALAME 
ASOACL*ATDET ( GAMMA , EM , ASTAR ) 

ASOATL*ASOACl + AL.OG ( AC / AT > 

C CHANGE ALAME TO DRIVE L Oiii ( AS TAR/AT ) TO ZERO 
CALL I TSUB ( ALAME , ASOATI. . . 20 , . 000 1 , SA VF 1 1. ) 

IBR»IFIX(SAVEIL( 1 ) ) 

GOTO < 20 , 20 , 20 , 20 , 300 , 30' » ) , I BR 
300 CONTINUE 

C SOLUTION HAS CON'TRGED -- PRINT COMBUSTOR AND THROAT RESULTS 
RH0*P*W/(C1*T) 

VEL*1. /<RHO*AC) 

WRITE (6, 9930) 

WRITE <6, 9800) 

WRITE(6,9810) T, P, RHO, FM, VFL 
WRITE(6,9840) W 
WRITE(6,9410) 

WRITE<6.9820) ( XA( I ) , 1=^1 , 6 ) 
li«ITE(6f9830) (EMUA(I),I=^1,6) 

WRITE<6»9100) . 

TS»T*TFOT/ ( . b* ( GAMMA+ 1 . ) ) 

PS»P# ( ( TS/T ) ** < GAMMA/ < GAMMA ~ 1 . ) ) > 

RHOS*PS*W/ ( C 1 *TS > 

VEL,S«1./<RH0S#AT> 

AISPS- (PS*AT+C3*VELS ) / ( 9 . 807*C3 ) 

WRITE<6.9900) 

WRITE <6, 9910) AISPS,TS,PS,RHOS,VELS 

WRITE 9100) 


ENTER THE CALCULATION OF THE EQUILIBRIUM EXPANSION TO THE THROAT 





o o 


C SAVE IDEAl SOI ur ION f I'f- INITIAI GUrSS 
1000 OAMMAC»GAMMA 
TC:«T 
PC*P 
FC-PF*AC 

RHOC=RHO 

INITIALIZE LOOP IN WHICH LMU(OH) WILL iJE VARIED TO 
CONSTRAINT 
SAVEIL(l)=l. 

SAVEIL(2)=-. 4*EMUA(5) 

L00P4«0 
DO 1010 1=1,6 
EMUAC( I )=EMUA< I ) 

1010 CONTINUE 
1020 DO 1030 1=1,6 

•EMUAl I )=EMUAC< I ) 

1030 CONTINUE 
ANM1=AC 
GAMMA^-GAMMAC 
FNMl-FC 
RHONMlsRHOC 
PNMl-P 
U«WC 
AN=AC 
T«TC 




SATISFY THE THROAT 


P=PC 

C DELTA MAY BE CHANGED TO MALE SOLUTION MORE OR LESS ACCURATE 

DELTA-. 01 
L00P3-0 

C REDUCE P FOR EXPANSION - NOTE THAT THIS IS IGNORED FOR FIRST PASS 
1100 P-P*<1. -DELTA) 

1110 IF <L00P3.£Q.0) GOTO 1120 
QA».5»(P-PNM1 ) 

QB= . 5* ( P+PNM 1 ) #ANM 1 -FNM J 
QC»C3*CUT/(P»W) 

C COMPUTE THE DUCT AREA THAT WILL SATISFY THE MOMENTUM EQUATION 

C KNOWING THE PRESSURE - IGNORE ON FIRST PASS 

AN-. 5* < -QB-SQRT ( QB*0B-4 . * OA«QC ) ) /QA 
GOTO 1130 
1120 QB-FNMl/AC 
QA-1. 

QC-C3»C1 *T/ < AC*AC*W ) 

C COMPUTE THE PRESSURE THAT WILL SATISFY MOMENTUM KNOWING AREA 
C FIRST PASS ONLY 

P« . 5# < QB+SQRT < QB*QB-4 . * OA^PC » ) 

C COMPUTE FORCING FUNCT I UH 

1130 BE-C4«C1*C1/(AN«AN) 

BM«C3«C1/(AN*AN) 

ARATii. 5« ( ANMl 4-AN ) / AN 
ENU-T/<P*W> 

BMENU-BN*ENU 
OP-SHENU/ < P »ARAT-BMENU ) 

PP10P^*ARAT/ < P#ARAT“ BMENU ) 


II 


oo 


C GET NEW THERMO PROPER nr:. 

CALL HASA 

XNUM 1 «BE<»ENU*ENU*Ff • 1 Of- -C OP/ W 

XNUM2=T*DHDT+BE^rNlJ*FI'IH»PP toP 
ALAME“-XNUM 1 / XNUM2 
DO 1140 1 = 1,/^. 

PHI ( I )»ALAME*(HA< I ) /T-Wi<^|iHD7 ) 
1140 CONTINUE 




ICNST«5 

GET NEW EQUILIBRIUM SOl.UTION FOR THIS nilCT POSITION WHILE HOLDING 
EMUAIOH) CONSTANT 
CALL EQUILdCNST) 

CALL HASA 
ENU»T/(P*W) 

DHFDT=-DHDT-BE f ENU*ENU/ I 
DHF«HF-H-. 5*BE*ENU*ENU 
T=T-.5*DHF/DHFDT 
C DRIVE ENERGY ERROR TO /( Rn 

IF <ABS(DHF).GT. DELTA) GOTO 11 tO 

C COMBUSTOR SOL UTION FOR THIS Mi:iL F FRACTION VAIUE HAS CONVERGED 
RH0»•P*W/<C1^^T) 

VEL=1./<RHCh«AN) 

DL0GP=2 . tt ( PNM 1 -P ) / ( PNM 1 +P ) 

DL0GR=2. * ( RHONMl -RHO ) / • RHONM I < RHO ) 

IF (L,00P3.GT.0) GAMMA HI. I it.l 

EM=SQRT ( BM*T/ ( GAMMA* W i ) /P 
I F ( LOOPS . OT . 0 ) GOTO 1 1 • 

C SAVE CHAMBER RESULTS 

SC=S 
TC«T 
• PC«P 
WC»W 

RHOC=RHO 

EMC«EM 

VELC-VEL 

DO 1150 1*1,6 

EMUAC(I)»EMUA<I) 

1150 CONTINUE 
1160 L00P3-L00P3+1 
RHONM 1-RHO 
ANMl-AN 
PNMl-P 

PNM 1 • < P4BM*ENU ) * AN 
A80ANL-ATDET ( GAMMA , EM , AC.OAN . 

C SINCE THROAT IS A SINGULARITY FOR THIS FORMULATION GET THE PROPER 
e AREA WITHIN 1 PER CENT AND USE IDEAL GAS EXTRAPOLATION BEYOND THAT 

C IF NOT WITHIN TOLERANCE GO BACK AND REDUCE PRESSURE 

IF <AB8{AS0ANL).GT..01) GOTO 1100 
A8TAR»A80AN»AN 
A0ATL**A80ANL-t>AL00 < AN/ AT ) 

C ALTER AMOUNT OF OH WHICH Will. BE PRESENT TO SATISFY THROAT 
CALL IT8UB(EHUAC<5>,A0ATL,20. .OOl.SAVEIL) 

IBR»lFIX<SAVEIL<t>> 

OOTO < 1020, 1020, 1020, 1020, 1200, 1200) , IBR 



C SOLUTION HAS COrWERGm . PfUNT RESULTS 
1200 WRITE(-i.,99?r)) 

WRITE<6,9800) 

WRITE(6,9810) TC, PC, RHi if , LML ,VFI.C 
URITE<A,984o) WC 
WRITE<6.9410) 

DO 1210 1-^1./. 

XA( I )=WC*FMUAi :( 1 ) 

1210 CONTINUE 

WRITE<6,98:-n) (XA<I),I-1 
WR I TE ( 6 , 98 .:iO ) ( EMUAC ( I ) , I 1 , ^ ) 

WRITE (6, 9 100) 

DO 1220 1=1.4 
XA(I)=W*EMUA(I ) 

1220 CONTINUE 

TS=T* < 1 . + . 5* ( GAMMA- 1 . ) *f. FH ' / ( . 5* ( Gi 'iMMA 1 1 . ) ) 

*PS=P»( <TS/T )•»••*( GAMMA/ (GAMMA- 1 . ) ) ) 

RHOS=PS*W/ ( C 1 *TS ) 

VELS«1./(RH0S«AT) 

A I SPS« < PS K AT4 C3^^ VELS )/ (‘K GOV F ?< ) 

WRITE (6, 9900) 

WR I TE < 4 , 99 1 o ) A I SPS , TS , P“: , RHGS , VEl. S 
WRITE (6, 9840) W 
WRITE<6,9410) 

WRITE(6,9820) (XA( 1 ), 1 = 1. ,4) 

WRITE(6,9830) <EMUA( I ) , 1=1 ,6) 

• WRITE(4,9100) 

STOP 

9 1 00 FORMAT ( 1 X , 1 30 < 1 Hfr ) ) 

9200 FORMAT< IX ,48HTHE RFSUl I' nf mi AnrAFtATfi M AMF ANALYSIS ARES) 

9300 FORINT ( IX. 14HTHE PRESSUlvE IS ,Gl3.6,i!X, 

* 19HTHE TEMPERATURE IS ,013.4) 

9400 FORMAT! IX, 23HTHE MOLE FRALflUNS ARE , 6( G13. 4, 2X ) ) 

9410 FORMAT! IX, 28X,2HH2, 13X,3HH20, 12X, IHH, 14X, IHO, 14X,2H0H, 13X,2H02) 
9800 FORMAT! IX, 35HTHE POST COMDUiniON CONDITIONS ARES) 

9810 FORMAT! IHO, 14HTEMPERATURE , G1 3. 4, 3X , 1 IHPRESSURE = ,G13.4,3X, 

* 4HRH0 = ,013.4, ,iiX, IIHMACH NO. = .G13.4.3X. 

* IIHVELOCITY = ,G13.4) 

9820 FORMAT! IX, 23HTHE M0L6 FRACTIONS ARE ,4(G13.6,2X) ) 

9830 FORMAT! IX, 24HTHE NUMBER OF MOLES ARE ,6<G13.6,2X) ) 

9840 FORMAT! IX, 24HTHE MOLECULAR WEIGHT IS ,G13.6) 

9900 FORMAT! IX, 26HTHE THROAT CONDITIONS AREs ) 

9910 FORMAT! IHO, 10HH1PULSE = , G1 3. 4, TX , l4HTEMPr.»!ATURE « ,013. 6, 3X, 

* 1 IHPRESSURE = ,G13.6,3X,4HRHG » ,013.6,3X, 

* IIHVELOCITY = ,013.6) 

9920 FORMAT! 1H1,51HTHE RESULTS FOR EQUILIBRIUM EXPANSION TO THE THROAT) 
9930 FORMAT! 1H1,49HTHE RESULlL. FOR IDEAL GAS EXPANSION TO THE THROAT) 
END 

FUNCTION ATDET!OAMMA, EM ASOAN) 

0M102*.5#!0AMMA”l. ) 

0F1-.5*(0AMMA+1 .)/ (GAMMA- 1 . ) 

0F2»AL00!.5«!0AMMA^1. ) ) 

VAL«AU00!EM)+GFHKGF2-AL0G< 1 . +GM102*EM*EM) ) 

ASOAN«EXP!VAL) 




c 

c 


ATDET=VAL 

RETURN 

END 

SUBROUT I NF F i M l 1 L < I C NS T ) 
DIMENSION ALAMJC) 




COMMON /SPf TiAT/f:MUA( . . I lA • Si ■ AA A ) , Al r, w . I , DHDT A < A ) , , 

* AJ(2) ,PH1 

COMMON /C0NGT/ATCU<2,6) < >u.:, /",/.) , WA(.;a , i. 1 , C2, C3 , C'I 
COMMON /STATE/T , P, W , H , : 

COMMON / M I SC / AC , AT , OAMMA 
LOOP 1=0 

10 ALAMJ(1)=-~.S*03A( 1 ) + PHI >; 1 )) 

ALAM J < 2 ) =- < GA ( 2 ) -GA ( 1 ) +PH I ( .• ; PH 1 U ) ) 

DO 20 1=3,6 

IFd.EQ. ICNS7 ) GOTO : > 

ARG=- < GPA ( I ) +PH I< I ) +AI. F.l-I- A J m>. ATi 0 U , I ) i /'.I AM J ( 2 ) * ATCO s 2,1 ) ) /C 2 
EMUA ( I ) = ( EXP ( ARG ) ) / < FaAJ > 

20 CONTINUE 
SUM«0. 

DO 30 1=3,6 

SUM=SUM+EMUA< I ) is ATCO <2, J ) 

30 CONTINUE 

EMUA < 2 ) = ( A J < 2 ) - SUM ) / AT L n ■: ) 

SUM=0 

DO 40 1=2,6 

SUM=SUM+EMUA ( I )*ATCD( t , I ' 

40 CONTINUE 

EMUA < 1 ) » ( AJ ( 1 ) -SUM ) / ATI. LA I , ll 

SUM»0. 

DO 50 1=1,6 

IF <EMUA( I ) . LT . . OOOCiAul > I MUa*. I )- . ai.ioOo'.U 
SUM=SUM+EMUA(I) 

50 CONTINUE 
W=1./SUM 
DO 60 1=1,6 
XA(I)»W#EMUA(I) 

60 CONTINUE 

LOOPl-LOOPl+1 

IF (LOOPl.LT.S) GOTO* 40 

RETURN 

END 

SUBROUT I NE HSOF T ( SPA ) 

DIMENSION SPA(6). 

COMMON /SPEDAT/EMUA ( 6 ) , tTA ( a ) , i lA ( 6 ) . GF A < /. ) , DHDT A ( 6 ) , XA ( 6 ) , 

* * AJ<2>,FH1<6) 

COMMON /CONST / ATCO < 2 , 6 ) , AA 2 . / , 6 ) , WA ( 6 ) , C 1 , C2 , C3 , C4 . C5 
COMMON /STATE/T , P , W , H , S , .P , iJl IDT 
COMMON /M I SC /AC , AT . GAMMA 

COMPUTE, GIVEN A TEMPERATURE, THE ENTHALPY AND ENTROPY <AT 1 ATM) 
FOR EACH OF THE SPECIES. 

KR»l 

IFIT.LT. 1000. ) KR=2 
DO too 1=1,6 
AI«A<KR,1, I) 


■ Mi 






S, Sj; - r 



c 


A2«A(KR.2. I ) 

A3«A(KR,3, I) 

A4-A(KR.4»I> 

AS»A(KR.S, I ) 

A6»A(KR.6,I) 

A7*A(KR»7, I) 

HP-Al-M A2/2. + ( A3/3. + ( A4/4 . + ( A5/5. ) *T ) *T ) *T > *T+A6/T 

SP«Al»AL0G<T)+(A2+(A3/2.+(A4/3.+(A5/4. )*T>»T)*T)*T 

SP«SP^A7 

HA<I)»C2*T*HP 

SPA(I)«C2*SP 




DHDTA ( I ) »C2* ( A 1 + < A2+ ( A3+ ( A4+A5*T ) *T ) *T ) #T ) 

100 CONTINUE 
RETURN 
END 

SUBROUTINE HAS A 
DIMENSION SPA(6),SA(6) 

COMMON /SPEDAT/EMUA(6) ,HA(6) ,GA<6).GPA<6),DHDTA(6)»XA<6) , 

* AJ(2),PHI'.6) 

COMMON /CONST/ ATCO( 2, 6 ) . A ( 2, 7, 6 ) , WA( 6 ) , Cl , C2» C3. C4. C5 
COMMON /STATE/ T . P , W . H, S , SP , DHDT 
COMMON /MISC/AC, AT. GAMMA 

ENTHALPY AND SYSTEM ENTROPY ARE DETERMINED KNOWING T.P.EMUA 
CALL HSOFT(SPA) 

H*0« 

S-0. 


SUM-0. 

DHDT-0. 

DO 100 1-1.6 
SUM-SUM+EML1A( I ) 
100 CONTINUE 
W-l./SUM 


DO 200 1-1.6 
H-H-»-EMUA<I)4HA(l) 

. SA ( I ) -SPA < I ) -C24AL0G < W*EMUA < I ) *P ) 

OA(I)-HA<I)/T-SA(I) 

OPA( I )-HA( I ) /T”SPA( 1 ) 

S-S-*>EMUA( 1 )4SA( I ) 

DHDT-DHDT-t>EMUA( I )*DHDTA< 1 ) 

200 CONTINUE 
RETURN 
END 

SUBROUTINE DATA^N 

COHNON /SPeDAT/EMUA(6).HA<6>.0A(6>.0PA(6)«DHDTA(6).XA(6). 
* AJ(2).PHI<6) 

COMNON /C0NST/ATC0<2»6).A(2.7.6).WA(6)>Cl»C2.C3.C4tC5 

COMMON /8TATE/T.P.W.H.S.SP.DHDT 

COMMON /MI8C/AC.AT.0AMMA 

READ(9t9IOO) PP.TF 

READt5*9200) <XA( I ) . 1-1.6) 

U-0. 

DO 100 1-1.6 

w-kH>xAm*WAa) 

100 CONTINUE 


I 

I 

} 

i 










t 


DO 200 1-1.6 
EMUA(I)-XA(I )/U 
200 CONTINUE 
T-TF 
P-PF 

CALL NASA 
HF-M 
SF— S 

WRITE(6.9300) HF.SF 
WRITE(6.9400) W 
DO 400 J-1.2 
AJ< J)-0. 

DO 300 1-1.6 

AU( J)-AJ< J)4-ATC0( J. 1 )*EMUA< I > 
300 CONTINUE 
400 CONTINUE 

MRITE<6.9500) ( AJ( J) , J-1 . 2) 




READ(S.9600> AC. AT. GAMMA 
WRITE<6.9700) AC. AT. GAMMA 
RETURN 

9100 F0RMAT<2E10.2) 

9200 F0RMAT<6E10.2) 

9300 F0RMAT<IH0.22HTHE TOTAL ENTHALPY IS .G13.6.3X. 

* 21HTHE FLAME ENTROPY IS ,013.6) 

9400 FORMAT ( IX, 24HTHE MOLECULAR WEIGHT IS .013.6) 

9900 FORMAT ( IX. 25HTHE NUMBER OF H ATOMS IS ,G13.6,3X. 

* 25HTHE NUMBER OF 0 ATOMS IS .013.6) 

9600 FORMAT (3E10. 2) 

9700 FORMAT! IX. 22HTHE COMBUSTOR AREA IS ,G13.6,3X. 

* 19HTHE THROAT AREA IS .G13.6.3X. 

* 27HTHE I8ENTR0PIC EXPONENT IS ,013.6) 

END 

SUBROUTINE DAT8ET 

. COMMON /SPEDAT/EMUA(6).HA(6).0A(6).0PA(6).DHDTA<6).XA(6). 
« AJ(2).PH1(6) 

COMMON /C0NST/ATC0<2.6),A(2.7,6).WA<6).C1,C2.C3.C4.C5 
COMMON /STATE/T.P.W.H.S.SP.DHDT 
COMMON /MISC/AC. AT, GAMMA 


C 

c 

c 

c 

t 


rd- ■ 



V 


SET THE ATOMIC COEFFICIENT^ FOP THE MOLECULAR SPECIES 

THE ARRAY ATCO<J.l) CONTAINS TH' NUMBER OF ATOMS OF ATOMIC SPECIES U 

CONTAINED IN MOLECULAR SPECIES I. THERE ARE SIX MOLECULAR SPECIES 

TREATED. THEY Af^E (IN ORDER) H2.H20.H.0.0H.02. ATOMIC SPECIES ARE H,0. 

ATC0(i.3)-l. 

ATC0(2f3)-0. 

ATC0(i,l)-2. 

ATC0(2fl)-0. 

ATC0(l*2)-2. 

ATC0(2.2)-I. 

ATC0a«4)«0. 

ATC0<2«4)«t. 

ATC0(it9)-l. 

ATC0(2t9)-I. 

ATC0(IiA)»0. 
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ATC0<2»6)=2. 

C IN THE FOLLOWING BLOCK THE THERMODYNAMIC COEFFICIENTS FOR EACH OF THE 
C MOLECULAR SPECIES ARE SET. THEY ARE STORED IN A(K,J.I) WHERE J«l.7 FOR 

C EACH OF THE MOLECULAR SPECIES 1=1,6. THE K INDEX IDENTIFIES THE 

C TEMPERATURE RANGE ( l=HIOH, 2-=LOW) . 

C H DATA FOLLOWS 

A(2, l,3)».25E-f01 
A(2,2,3)»0. 

A(2,3,3)=0. 

A(2,4,3)-0. 

A(2.S,3)=0. 

A<2,6,3)».2547l627E+05 

A(2,7,3>«=-.46011762 

A<l,l,3)=.25E+01 

A(l,2,3)=0. 

A( 1,3,3) =0. 

*A(1,4,3)=0. 

A(l,5,3)-0. 

A ( 1 , 6 , 3 ) = . 2547 1 627E-^05 
Ad, 7, 3)=-. 4601 1763 
C H2 DATA FOLLOWS 

A ( 2 , 1 , 1 ) » . 3057445 1 E-^0 1 
A(2,2, l)=.267652E-02 
A < 2 , 3 , 1 ) =- . 58099 1 62E-05 
A<2,4, l)=.55210391E-08 
A<2,5, 1>«“. 18122739E-11 
A ( 2 , 6 , 1 ) — . 98890474E+03 
A ( 2 , 7 , 1 ) «“ . 22997056E-^0 1 
Ad, I, I )».31001901E-^01 
Ad,2, l)-.51119464E-03 
Ad , 3, 1 ) ». 5264421E-07 
A d , 4 , 1 ) «- . 34909973E- 1 0 
A (1 , 5, 1 ) - . 36945345E- 1 4 
Ad,6, 1>— .87738042E+03 
Ad,7,l)»- .19629421E+01 
C H20 DATA FOLLOWS 

A<2, 1,2)«.4076i 275E4-01 
A(2,2»2)«-. 11084499E-02 
A(2,3,2)-.41521l80E-^5 
A< 2, 4, 2 ) *•*. 29637404E-08 
A(2,5»2>-.80702103E“12 
A<2,6,2)»-.30279722E+05 
A<2.7,2)*-.3227p046 
Adf i,2)«.27167633E-^01 
Ad f 2, 2)-. 29451374E-02 
A< 1 f 3, 2)*-. 80224374E~06 
A<1.4f 2)«. 10226682E>09 
Adt8«2>— .48472145E-14 
Ad * 2) — . 29905826E-*>05 

A d * 7 * 2 ) - . 6630567 1 E't’O I 
e 0 DATA FOLLOWS 

A(2, lt4>*.29464287E4^01 
A(2i2»4)— . t6381665E-02 
A<2.3*4>«.242103l6E-05 



ooo ooon 




RETURN 

END 

SUBROUTINE ITSUB(X,FOFX,NITSUB,CONV,SAVEIT> • 

DIMENSION SAVEIT(8) 

ITERATION CONTROL SUBROUTINE ...FINDS X SUCH THAT F(X)-0 
MAIN PR(M3RAM MUST DIMENSION SAVEIT<8) 

CALLING PROGRAM COMPUTES < GIVEN AN X) FOFX AND ALSO 
SUPPLIES CONVERGENCE CRITERION (CONV) AND « OF ITERATIONS 
(NITSUB)...ITSUB THEN CONTROLS SUBSEQUENT VALUES OF X 
INITALIZE ITSUB BY SETTING SAV^IT(1)«1 
SET AMOUNT TO PERTURB X BY IN SAVE I T< 2) 
IF((AreCFOFX>-CONV>.LE.O. ) GOTO 730 
ITIME«IFIX(SAVEIT(1)4-.1) 

GOTO (330,380,470,570). ITIME 
330 CONTINUE 
1TIME*2 
*SAVEIT(3)»I. 

FOFXCK-FOFX 
SAVEIT(8)-F0FXCK 
IF(FOFX.LT.O. ) GOTO 470 
3B0 CONTINUE 

IF(FOFX.LT.O. ) GOTO 570 
IF(FOFXCK.OE.FOFX) GOTO 430 
8AVEIT<2)»- SAVEIT(2) 

X-X-2.*SAVEIT(2) 

GOTO 670 
430 CONTINUE 
8AVEXT<4)»X 
8AVE1T(5>-F0FX 
X-X-8AVEIT<2) 

GOTO 670 
470 CONTINUE 
ITXME-3 

IF(FOFX.OT.O. ) GOTO 570 
IF(FOFXCK.LE.FOFX) GOTO 530 
8AVEXT(2>— SAVEIT(2) 

X«X^2.«8AVEIT(2> 

GOTO 670 
830 CONTINUE 
8AVB1T(6»X 
, 8AVEXT(7>«F0FX 
X«X43RVEXT(2> 

GOTO 670 
870 CONTINUE 
ITIMEM 

IFCrarX.LT.O. ) GOTO 630 
(HMCIT(4)«X 
8RVEXT(8><i>F0FX 
GOTO MX) 

630 OONTilNUE 
•RWEIT(6)«X 
•RVilT(7)-F0FX 
683 CONTINUE 

X«3RWBlT<4>-3RUElT<8>ft(8AVEtT(6>-6AVEIT(4> )/(8AVEIT(7>-8AVEIT(S> ) 


^ Won «TV 


1F(<X-8AVE1T<4>)*(X-SAVEIT<6)).0T.0. ) GOTO 710 

670 CONTINUE 

NTEHP«IFIX<SAVEIT<3>+. 1 ) 

IF(NTEMP.OE.NITSUB) GOTO 710 

NTEMP-NTEHP-H 

8AVEIT ( 3 ) -FLOAT ( NTEHP ) 

OOTO 780 
710 CONTINUE 
ITIHE-6 
OOTO 780 
730 CONTINUE 
ITIME-9 
8AVE1T(4)>X 
SAVEIT(5)-F0FX 
8AVEIT<6)-X 
8AVEIT(7)-F0FX 
780 t^ONTXNUE 

8AVEIT< 1 )-FLOAT( ITIME) 

RETURN 

END 






